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4.1 Introduction 

The first modern generation of signal processing algorithms used Fourier 
methods. Their theoretical basis was provided by the convolutions which 
occur in time-invariant linear systems and the important role of sinusoids 
as eigenfunctions of every such system. A second class of signal analysis 
and parameter estimation methods use a general linear model similar to the 
one used in the statistical literature. Representative methods in this class 
are the linear predictive and signal subspace methods for parameter esti- 
mation. A third class of signal processing methods based on time-frequency 
analysis of nonstationary signals began with the work of Wigner (1932) and 
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Ville (1948) and has seen extensive theoretical development and the begin- 
ning of hardware implementation in the 1978-1988 period. Recent work 
has combined linear predictive modelling with YVigner-Ville analysis. The 
near future wiil see the combination of time-frequency distribution (TFD) 
methods with both eigensystem-based signal analysis and wideband ambi- 
guity function and wideband WVD analysis in which Doppler is treated as 
a time compression or expansion. 

This paper provides an updated discussion of high resolution process- 
ing techniques for temporal and spatial signal processing, with particular 
emphasis on methods applicable to underwater acoustic signal processing. 
The paper consists of three parts: The first part treats direction finding 
and spectral estimation for stationary random processes. The second part 
discusses time-frequency analysis for nonstationary signals using both the 
Wigner-Ville distribution (WVD) and an autoregressive modified WVD. 
The third part discusses the application of these techniques to radar and 
sonar imaging. 



4.2 Direction Finding 

This section provides an overview of selected beamforming and direction 
finding (df) techniques and their computational requirements. It provides a 
description of algorithms and an assessment of their range of applicability. 
Last it recommends areas for future study. With a few noted exceptions, 
it considers only narrowband processing - i.c. it is assumed that the data 
snapshots are the Fourier coefficients in one frequency bin for each of the 
array elements for the currrent observation interval, and that different fre- 
quency bins are processed independently. 

4.2.1 Algorithms Considered 

• Classical Beamforming 

• MVDR and Capon’s Method 

• DF by Linear Predictive Spectral Estimation 

• MUSIC 

• Johnson & DcGraf 

• Other MUSIC variants 



* Pisarenko Harmonic Retrieval 
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• ESPRIT 

• ESPRIT variants 

• Cadzow’s Eigenvector Approximation Method 

• Maximum Likelihood 

4.2.2 Principal Notation for Direction Finding 

x data vector 

data vector for kth snapshot 
Iv total number of snapshots currently available 

M number of elements in array 

B number of beams 

W total processed bandwidth 

a(u) arrival phase shift vector for look direction 11 . 

As a function of u, also called the “array manifold” 

A 1) general matrix 

2) steering phase shift matrix A = (a(ui) • • ■ a(uo)) 

D 1) number of arriving wavefronts (can be larger than 

the number of sources) 
when multipath arrivals are considered) 

2) diagonal matrix in svd or eigenvalue decomposition 
.s steering vector, alternate notation for a(u) 

E usually denotes statistical expectation 

II used as a superscript to denote Ilermitian transpose 

P 1) source correlation matrix 

2) unitary matrix 
Q 1) unitary matrix 

2) matrix with orthonormal columns - part of a unitary matrix 
U usually used to denote upper (or right) triangular matrix. 

In the numerical analysis literature, R is usually used to denote 
a right triangular matrix. However, since the signal processing 
literature usually uses R for a correlation matrix, we will usually 
use U instead, except in the names of numerical algorithms, 
such as QIl decomposition and the QR cigcnsystcm method. 

R spatial correlation matrix. R = Exx" 

When it is desired to emphasize the dependence on x, we will 
use /7.x 

R n 1) noise spatial correlation matrix 

2) scaled noise correlation matrix 
w weight vector (combined ampii tude and phase) 
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PDQ fI factors in singular %’aluc decomposition. P is unitary 
D is diagonal, and Q is unitary 

QDQ 11 factors in eigenvalue decomposition of a Ilermitian matrix. 

Q is unitary, D is diagonal, 
n noise vector 

f 1) vector of complex amplitudes for wavefronts arriving at 

the array 
2) frequency 

<r 2 noise power, possibly after spatially prewhitening 



4.3 Overview and Algorithm Assessment 

4.3.1 Classical Beamforming 

We consider classical beamforming only to establish notation and termi- 
nology that will be used in discussing additional algorithms. Standard 
references for classical beamforming include Ma [30] and Steinberg [65], 
For most algorithms, we will regard the frequency-domain output of the 
current data snapshot for the array as a single complex vector, x. The data 
model used is: 

x = /If + n 

Classical frequency domain beamforming via phase shift and sum is the 
frequency domain equivalent of delay and sum time domain beamforming. 
For each look direction u, it forms the inner product of the steering phase 
shift vector s = a(u) with the data vector x. This inner product, b = 
[a(u),x] is the spatial equivalent of a matched filter. A derivation for general 
array geometries with simplifications for special array geometries has been 
provided by Speiser, Whitehouse, and Berg [59]. The function a(u) is 
sometimes called the “array manifold” [52]. 

Shading may be used to reduce sidelobe height, at the expense of re- 
duced resolution. Special array geometries permit the beamforming to be 
performed via an FFT. These include uniform line and planar arrays, cir- 
cular and cylindrical arrays, and logarithmic linear and spiral arrays [59]. 

4.3.2 MVDR Beamforming 

Minimum Variance Distortionless Response (MVDR) beamforming attempts 
to reject interfering sources while maintaining unit gain and zero phase shift 
for each look direction MO] [43]. It is sometimes referred to a,s “Optimum 
Array Processing” [40]. Unlike classical shaded beamforming, it does not 




133 



try to create uniformly low sidelobes in the response pattern, but only to 
form nulls in the direction of interfering sources. For each steering vector 
s = a(u), it forms b = [w,x], where the weight vector w is chosen to mini- 
mize the expected power output, while providing unit gain and zero phase 
shift for a signal arriving from direction u. i.c. pick w to minimize E\b\ 2 
subject to [w,s] = 1. 

R-'s 

Wopi ~ [R~ l s t s] 

MVDR is not a high resolution method per se, but it is recommended 
as a baseline for comparison. It is quite robust, useful with any geometry, 
and very well understood, both in terms of its performance and parallel al- 
gorithms and architectures for real time implementation, including variants 
incorporating multiple linear constraints [54] [32], 

Apart from resolution MVDR does have one limitation: coherent mul- 
tipath arrivals in the same frequency bin may cause it to effectively lose 
sensitivity. That is, when the beam is steered to one of the arrivals, even 
though the weight vector provides unit gain in that direction, the response 
to a coherent multipath arrival may reduce the signal component of the 
output. 

4.3.3 Capon’s Method 

Capon’s method estimates the angular spectrum, or power versus arrival 
direction as l/(/2 _1 s,s] [17], [18], [26]. This is the average power output 
of an MVDR beamformer with steering vector s. The signal processing 
literature (including Schmidt’s papers) frequently refers t.o Capon’s method 
as “maximum likelihood”. This is a source of confusion, since it does not 
perform maximum likelihood parameter estimation, and does not have the 
high resolution of true maximum likelihood spectrum estimation/direction 
finding. We will regard Capon Maximum Likelihood as just an auxiliary 
output easily available from an MVDR beamformer. 

4.3 A DF by Linear Predictive Spectral Estimation 

Linear predictive spectral estimation determines the weights of a prediction 
error filter for a wide-sense stationary random sequence. 

Such a filter converts the sequence into white noise. Since the spectral 
density of the output of a filter is the product of the input, spectral den- 
sity with the magnitude squared of the filter’s transfer function, the input 
spectral density function is proportional to the reciprocal of the magnitude 
squared of the filter’s transfer function. When the input process has an 
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all-pole spectrum, the optimum prediction based on the infinite past is the 
same as the optimum predictor using a number of past samples equal to 
the number of poles. This method is sometimes called Maximum Entropy 
Spectral Estimation - more specifically the Yule- Walker version of Max- 
imum Entropy Spectral Estimation [26], [35], It has been proposed as a 
direction finding technique by utilizing the equivalence of spectral estima- 
tion for a wide-sense stationary random process with power versus arrival 
direction estimation for a uniformly spaced line array receiving signals from 
uncorrelated sources [31]. Linear predictive spectral estimation can provide 
high resolution, but require special precautions to avoid line splitting and 
bias errors [35]. The optimization of a signal processing function analogous 
to a thermodynamic function can be used not only to derive linear predic- 
tive estimators, but also can be used in a new class of estimators which can 
be used with a small number of data snapshots [56]. 

4.3.5 MUSIC 

The MUSIC method is at this time the most thoroughly studied and best 
understood of the eigensystem based direction finding techniques. It was 
independently developed by R. Schmidt [52], [53] in the U. S. and by the 
team of L. Ivopp and G. Bienvenu in France [6]. Schmidt’s formulation 
uses the generalized eigensystem of the data correlation matrix R x and the 
scaled noise alone spatial correlation matrix 7?n ■ The Bienvenu and Kopp 
formulation (and many others who have studied MUSIC) assume that the 
noise is spatially white or the data has been spatially prewhitened, so that 
Rn = o -2 /. The noise power, <r 2 is regarded as an unknown parameter. In 
this case, the generalized eigensystem reduces to the ordinary eigensystem 
of R x alone. In discussing MUSIC it is necessary to distinguish between the 
behavior with perfect knowledge of the data correlation matrix and scaled 
noise alone correlation matrix, versus behavior with estimated or modelled 
correlation matrices. 

MUSIC requires several assumptions: 

1. The number of arriving wavefronts (in the frequency bin being ana- 
lyzed) is strictly less than the number of elements in the array. 

2. The columns of the arrival phase shift matrix, A, are linearly inde- 
pendent. 

3. The source correlation matrix, P, is strictly positive definite, (sources 
can be coherent - but not perfectly coherent) 

4. The noise spatial correlation matrix is nonsingular. (This requirement 
is unlikely to ever cause a problem) 
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With perfect knowledge of the data correlation matrix and scaled noise 
spatial correlation matrix, MUSIC proceeds as follows: 

Solve the generalized eigensystcm 



R x z = Af? n z 



If the generalized eigenvalues are numbered in decreasing order, the last 
M — D are equal, and the first D generalized eigenvalues are strictly greater. 
The span of the generalized eigenvectors corresponding to the last M — D 
generalized eigenvalues is called the noise subspace. Each vector in the 
noise subspace is orthogonal to the arrival phase shift vector a(u,i) from 
each source direction. The number of sources is estimated as D1 = M — Mo 
where Mo is the multiplicity of the smallest generalized eigenvalue. 

The arrival directions are then determined by looking for array mani- 
fold vectors a(u), orthogonal to the noise subspace. This determines the 
arrival phase shift matrix A. Then the source correlation matrix P, can be 
computed as 

P = {A lt A)~ l A" {R x - A mi nPn)A(A ,l A)- 1 . 

Of course it is preferable to use the pseudovinver.se of A , computed by an 
SVD, rather than forming 

{A" A)'* A n 

in order to avoid problems of ill-conditioning caused by near collinearity of 
columns of A when signals arrive from nearly the same direction. 

Early simulations by Schmidt, using perfect correlation matrices, showed 
that MUSIC had far higher resolution than classical beamforming and 
Capon’s maximum likelihood method, and did not have the bias problems 
of the maximum entropy method. 

It is important to note that the MUSIC “DOA spectrum” is simply a 
function that peaks in directions where the steering vector is nearly orthog- 
onal to the noise subspace. That is, peaks correspond to arrival directions, 
but the DOA spectrum is not a measure of arrival power versus direction. 
Furthermore, the processing is highly nonlinear, so there is no simple inter- 
pretation of the apparent pcak-to-sidelobe ratio. 

The previous comments must be modified somewhat when estimated 
correlation matrices are used, but MUSIC still remains a very strong con- 
tender for direction finding. When estimated correlation matrices arc used: 

1. Instead of a smallest generalized eigenvalue with multiplicity (M - 
D), there is a spread of small generalized eigenvalues. If wc pick too 
large a dimension for the noise subspace, then signals will lie missed. 
If we pick somewhat too small a dimension for the noise subspace, 
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the method degrades gracefully: Spurious peaks will be found in the 
DOA spectrum, but these can be eliminated through estimation of 
the source correlation matrix and eliminating peaks corresponding to 
zero power. 

2. With ideal data and scaled noise correlation matrices, MUSIC can 
localize in the presence of highly correlated (but not perfectly cor- 
related) multipaths. In simulations using estimated correlation ma- 
trices, problems have been reported when the correlation between 
multipath arrivals was around .7 or so. 

A survey of work on MUSIC and related eigenvector techniques through 
1985, with detailed discussion of source covariance matrix estimation and 
a signal estimation techniques is available [61]. 

4.3.6 Johnson DeGraaf 

The method of Johnson and DeGraaf is a slight modification of MUSIC. 
In MUSIC, the generalized eigenvalues are only used in deciding on the 
dimensionality of the noise subspace, not in forming the DOA spectrum - 
i.e. all the basis vectors for the noise subspace are weighted equally in the 
DOA spectrum. Johnson and DeGraaf use a modified DOA spectrum, in 
which the eigenvector weighting depends on the corresponding eigenvalue, 
so that in effect the decision as to the dimension of the noise subspace (or 
equivalently the number of signals present) is a soft decision [27]. 

In simulations by the inventors, it was found to provide the same perfor- 
mance as MUSIC when the correct number of signals was chosen for both, 
but to be more robust when errors were made in picking the number of sig- 
nals. Studies by Gordon Martin confirmed this behaviour when simulated 
data was used, but showed better performance for MUSIC when actual IIP 
electromagnetic data was used [36]. 

4.3.7 Other MUSIC variants 

Two MUSIC variants intended to permit use with highly correlated multi- 
path arrivals are applicable for uniformly spaced line arrays. One technique 
uses averaging of the correlation matrices from subarrays [55], [49]. In ef- 
fect this tends to restore the Toeplitz structure that one would expect with 
uncorrelated sources. However, it both reduces the effective aperture and 
the number of arrivals that can be handled by a factor equal to the num- 
ber of subarrays used. The second technique is called Root MUSIC. For 
a uniform line array, the steering vectors are sinusoids, so finding steering 
vectors orthogonal to the noise subspace is equivalent to finding zeros of 
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polynomials on the unit circle. In Root MUSIC, zeros are rejected if they 
lie too far from the unit circle, and zeros close to the unit circle are pro- 
jected onto the unit circle [50]. Root MUSIC is reported to provide better 
threshold performance than classical MUSIC. 

4.3.8 Pisarenko Harmonic Retrieval 

Two spectral estimation techniques are both sometimes referred to as ’’Pis- 
arenko’s Method”. One is a family of methods based on the eigenvalues 
and eigenvectors of the correlation matrix. Particular members of this 
family are chosen by picking a monotone function and its inverse. In this 
report, we discuss only the Pisarenko special harmonic retrieval method 
[48]. This method determines the frequencies and amplitudes of sinusoids 
in additive white noise. Although it preceded (and inspired) MUSIC, it is 
best understood as a special case of MUSIC, with the following additional 
assumptions: 

1. The noise is white. 

2. The data correlation matrix is Toeplitz. In df application this corre- 
sponds to uncorrelated sources and a uniformly spaced line array. 

3. The noise subspace is treated as having dimension equal to 1. 

The Pisarenko method requires finding the eigenvector or a Toeplitz cor- 
relation matrix corresponding to the smatlest eigenvalue, and then finding 
sinusoid vectors orthogonal to that eigenvector. Although it can provide 
high resolution it is not recommended for the following reasons: 

1. The Toeplitz assumption makes it inapplicable when even partially 
coherent multipaths are present. 

2. Simulation studies and theoretical analysis have shown that the fre- 
quency (or angle in the df application) estimates have large variance 
or poor statistical efficiency. 

3. Gordon Martin’s studies using IIP data have confirmed the large vari- 
ance of the direction estimates [36]. 

4.3.9 ESPRIT 

ESPR.IT requires an array consisting of two identical subarrays, with the 
same displacement vector between each element of the first array and the 
corresponding clement of the second array. Alternatively, the total set 
of array elements may be regarded as a set of M dipoles, with the same 
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differential displacement vector (length an orientation) for each dipole. The 
following advantages are claimed for ESPRIT in comparison with MUSIC 
[45], [46], [47], [50]: 

1. ESPRIT eliminates the computationally expensive search through the 
DO A spectrum. 

2. It provides better performance in the presence of highly correlated 
multipaths. 

3. It provides a lower noise threshold. 

4. It does not require knowledge of the array manifold. 

However, the following disadvantages of ESPRIT for spatial signal pro- 
cessing should also be noted [63]: 

1. It is applicable only to a very limited set of array geometries. 

2. Even when a two-dimensional array is used, it only produces the same 
information as other methods using a line array - i.e. arrival angle 
with respect to the differential displacement axis. 

3. If the length of the differential displacement vector exceeds half a 
wavelength, then ESPRIT will have grating lobe ambiguities. 

For applications other than spectrum analysis and direction finding with a 
uniformly spaced line array spaced closer than spatial Nyquist, the above 
problems limit the usefulness of ESPRIT. 

4.3.10 ESPRIT Variants 

The original formulation of ESPRIT was based on solving a non-symmctric 
generalized eigensystem. Two recent alternatives have been proposed: To- 
tal Least Squares ESPRIT, and Orthogonal Procrustes ESPRIT [78], [7 9], [80], 
[51]. Total Least Squares ESPRIT is claimed to provide slightly better per- 
formance at low signal to noise ratio. 

4.3.11 Cadzow’s Eigenvector fitting method 

Cadzow’s method is a modification of MUSIC that avoids the requirement 
that the source correlation matrix be strictly positive definite [16]. That 
is, perfectly coherent sources are permitted. The method starts with the 
eigensystem and “DOA spectrum” calculations as in MUSIC. The noise 
is assumed to be spatially white, or prcvvhiteiied if necessary. Then it is 
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observed that each eigenvector of the data correlation matrix is a linear 
combination of steering vectors from the source directions. Source direc- 
tions are then found by fitting each eigenvector by linear combinations 
of steering vectors. Usually only a few steering vectors (corresponding to 
highly correlated arrivals) are used to fit each eigenvector. Since the steer- 
ing vectors depend nonlinearly on the arrival directions, this is a nonlinear 
least squares problem, and must be solved iteratively. 

4.3.12 Maximum Likelihood Parameter Estimation 

This technique (not to be confused with Capon’s so-called “maximum like- 
lihood” method) is based on the statistical maximum likelihood parameter 
estimation method, and assumes that the noise distribution is Gaussian. 
Maximizing the likelihood ratio for jointly estimating the signal amplitudes 
and arrival directions is then equivalent to a nonlinear least squares fit of 
the data vectors by linear combinations of steering vectors - a computa- 
tional task similar to that required by Cadzow’s method. It is known from 
the statistical literature that maximum likelihood parameter estimation has 
excellent large sample properties - asymptotically unbiased, consistent, and 
asymptotically efficient. Recent empirical studies of maximum likelihood 
direction finding suggest that the small sample behavior is very competi- 
tive with other high resolution techniques, and that the method is robust 
in the presence of highly coherent sources [14], [72], [77], It also extends in 
a straightforward way to broadband processing. 



4.4 Introduction to Nonstationary Signal 
Analysis using the Wigner-Ville 
Distribution 

An increasing set of signal processing tasks require an extension of spectral 
analysis to treat nonstationary functions of a time or space variable. Tradi- 
tional techniques for the analysis of signals and stationary random processes 
are based on the fact that the sinusoids are eigenfunctions of every shift- 
invariant linear operator. This includes both the impulse response of a 
time-invariant linear filter and the covariance funct ion of a widc-scnse sta- 
tionary random process. The Karhuncn-Locvc (K-L) expansion represents 
an arbitrary random process in terms of the eigenfunctions of its covari- 
ance function, and therefore permits an expansion in terms of uncorrelated 
random variables. However, the K-L expansion docs not posess convenient 
properties for the analysis of signals that have been filtered or modulated. 
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Recent signal analysis techniques using time-frequency distribution (TFD) 
functions provide a bridge between the the traditional stationary signal 
analysis methods and the yet to be developed methods for treating general 
nonstalionary problems. The Cohen [19] class of time-frequency distribu- 
tions provide a time-varying spectrum with many useful properties for a 
wide class of signals. Different TFDs in the Cohen class are selected via 
the choice of a weight function. The choice of a weight function determines 
which useful properties will be posessed by the resulting TFD, and no one 
choice provides all of the possible desirable properties. 



4.4.1 Principal Notation for Time-Frequency Analysis 



t 

f 

«( 0 
n 

s(t) 

A" 

z(l) 

S*(LJ) 

*/ 

*t 

*</ 

m(<) 

fo 

fc 

h 

a 

T 

e(n) 

c(», j) 

«(0 

C+ 

A 



time 

frequency 

imaginary unit (also used as an index) 
real signal 
Ililbert transform 

alternate notation for the Ililbert transform of s(t) 

Hermitian transpose of a matrix A 

complex analytic signal 

Wigner-Ville distribution of z(t) 

short-time Fourier transform of z(L) 

bilinear kernel corresponding fo z(l) 

convolution with respect to frequency 

convolution with respect to time 

two-dimensional convolution in (,/ 

window or modulation function 

initial frequency 

carrier frequency 

sampling frequency 

chirp parameter 

length of analysis window 

white noise sequence 

estimated covariance matrix 

autoregressive model parameters 

Moore-Penrose pseudoinverse of matrix C 

eigenvalue 

eigenvector 

standard deviation (also used for singular value] 
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4.4.2 The Wigner-Ville Distribution 

One TFD that has a large set of the possible useful properties-, and that has 
received extensive study in recent years for application to signal processing 
is the Wigner-Ville Distribution (WVD) (73],[37). The WVD of a reai signal 
s(t) that has complex analytic signal r(£) is denned as 

W z (iJ) = J z(l + r/2)z*(t - T/2)e~ i2 *' T dr 

where z(i) — s(t ) + is(t) and s(t ) denotes the Hilbert transform of s(£). 
The WVD can also be expressed as 

W t (l,f) = J I< t (t,T)c- i2 *' T dr 

where the kernel K t (t, r) is defined as 

A\(f,r) = z(t -f r/?.)z'(l - r/2). 

The WVD of a chirp (a function with constant amplitude and quadratic 
phase) is distributed on a line through the origin of the time-frequency 
plane, with slope proportional to the chirp rate. The WVD of the convo- 
lution of two signals is the convolution in time of their respective WVDs: 

HW*,/) = J W z (t-uJ)W h {u,f) du 

These two properties are crucial to a number of applications of the WVD: 

• Location of Sources in the Fresnel Region 

• Coherent Integration of Lines Emitted by a Moving Source 

• Identification of Stationarity Intervals for a Nonstationary Random 
Process 

• Radar Imaging 

Since a narrowband signal or signal component in the Fresnel region of a 
line array produces a quadratic phase shift as a function of position along 
the array, such a source can be localized in range and bearing via a spatial 
WVD in the temporal Fourier transform domain. [13]. 

Estimation techniques using signal subspace concepts can be applied 
to time-frequency distributions by replacing the Fourier transform of the 
bilinear kernel by a parametric modelling with respect to the variable r at 
each time i [8], [12], Such a parametric high-resolution modified WVD can 
be used to enhance the resolution attainable in the spatial localization of a 
source in the Fresnel region of a line array [74], [G6]. 
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4.5 The Windowed Wigner-Ville Distribu- 
tion 

Let s(t) be a real continuous-time signal. Tlie analytic signal associated 
with s(t) may be expressed as 

z(t) = s(t) + i7i[s(l)} (4.1) 

Where Ti denotes the Hilbert Transform. 

The WVD of a signal s(t) is defined as follows: 

/ + °° 

z(< + r/2)z*(f-r/2)c"^ T (4.2) 

■ OQ 

where z(t) represents tlie analytic signal defined in (1). The importance 
of tlie use of this analytic signal is described in [11]. 

The properties of the WVD are well described in several comprehensive 
papers, such as [37] and [10]. In particular, the WVD exhibits a range of 
peaks along the curve describing the instantaneous frequency law of the 
signal, which allows frequency tracking of time-varying signals. 

Inspection of Equation 4.2 indicates that the evaluation of the WVD is a 
non-causal operation. Thus, the value of the signal should be known for all 
time before either the Hilbert transform or the WVD can be evaluated. In 
many practical cases, the signals considered have either a finite duration or 
band limited spectrum and the problem is overcome by applying the WVD 
to a windowed version of the signal z(l) at time £ 0 such that : m (l,t o) = 
z(£)m(£- to), where m(t) is a time-limited window function with duration T. 
The WVD of the windowed signal is given by a one-dimensional convolution 
in the frequency variable, /: 

W 2m (tJ)=WAtJ)* f W m (tJ) (4.3) 

where */ denotes frequency convolution. 

Thus the effect of windowing is to smear the WVD representation in the 
frequency direction only; hence the frequency resolution can be increased 
by using a longer window without adversely affecting the time resolution. 
This is a significant improvement over the magnitude squared short-time 
Fourier transform (STFT) S Zm (l,f) where time resolution is degraded by 
the window as shown by the following relationship 

\S,Jl,f)\ 3 = w m (I.J) 

where denotes convolution in time and frequency. 



(4.4) 
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For notational convenience, we define K z and k m such that: 



K*(t t r ) = z(t + r/2)z‘(l - r/2) and h m (r) = m(r/2)m* (-r/2) (4.5) 

The WVD at any point in time is evaluated by shifting the signal z(t) so 
that the window is centered about t — 0. Therefore, the WVD computation 
reduces to the evaluation of Equation 4.6: 

,+T/2 

W t Jt,f)= / K : (t,r)lc m (r)e- i2 * lT dT (4.6) 

J-T/2 

where W z is an estimator of W z . 

Equation (6) shows that the estimation of the WVD reduces to the com- 
putation of the FT of a signal kernel K 2 (1,t) as a function of r with an 
analysis window h(i) = lc, n (T), which is a classical problem of signal anal- 
ysis. An efficient implementation of (6) which uses fundamental symmetry 
properties of the WVD is described in [9]. 

4.6 The Autoregressive WVD (AWVD) 

4.6.1 Description of the method 

Equation (6) shows that WVD analysis may be thought, of as a two-step 
process: 

1. A sequence of kernels t) is formed. 

2. A Fourier Transform of /\%(<,r) is taken, with respect to the variable 
r, using the analysis window k m to reduce spectral leakage. 

It is easily shown that for the linear FM signal, A' 2 (f,r) becomes a sin- 
gle complex sinusoid with frequency varying with the chirp parameter tv 
determining the slope [7]. That is, if 

z(l) = e ‘2*(/cl+crl a /2)^ (4.7) 

then the WVD kernel becomes: 

Kg{t, t) = e i2r Uc+a()T (4.8) 

Since a harmonic process can be generated by an autoregressive (AR) 
model, the real part of the kernel r) at time t may be fitted by the 

output of an appropriate AR model; the AR coefficients will estimate the 
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cross-section of the Wigner-Villc Distribution at time t. The real part is 
taken because the kernel has an llermitian symmetry and therefore, no ad- 
ditional information is contained in the imaginary part. The analytic signal 
must however be used in the estimation of the kernel. The expression for 
the kernel is then: 



K.{L,t) = s(t + r/2)s*(t - r/2) + 7i[s(t + r/2)]?f[s*(t - r/2)] (4.9) 

where W(.) represents the Hilbert Transform. A variety of methods for 
super-resolution spectral estimation are available for use with this method, 
based on parametric modelling. 

i) Autoregressive modelling 

A method of spectral estimation which has been widely accepted in seismic 
[15] and speech processing [33] is the method of auto-regressive modelling. 
Makhoul [33] shows that the autoregressive modelling is equivalent to solv- 
ing the normal equations associated with the linear prediction problem for 
a finite set of observed samples. For the optimum predictor to use only a 
finite number of samples of the past, the prediction error filter will have 
only zeros, and the spectrum being estimated must have only poles. It is 
the use of this structural information which permits the use of maximum 
entropy estimators when the process model is appropriate [33]. 

The traditional approach for implementing this method consists of 

1. calculating the sample autocorrelation or autocovariancc function 

2. solving a linear least-squares problem to find the optimum predictor 
tap weights 

3. Fourier transforming the prediction error filter weights. 

This method assumes that the sampled kernel data under analysis can be 
expressed as a p th order AR model. Then the N point truncated kernel 
must satisfy the difference equation 



r 

M«) = f(u) - a (0M n ~ *') 

i=i 



where N is odd and 



M»») 



z{i -f- Ij)z (/. — Z/) 

n - r 



(4.10) 



(4.11) 



L 



n — 
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c(k) is a sample of the noise which is assumed to come from a white noise 
process; a(i), i — 1 ...p are the coefficients of the AR model. (Here, we 
have assumed a change of variable n = r/2 before discretizing, as reported 
in [11].) 

Then the covariance AR parameter estimates, a and &, arc given by the 
solution of 

Ca = <t 2 ci (T.l 2) 

where C is the covariance matrix formed by the elements 

1 JV_1 

c(»,i) = j) fori,j=0...p (4.13) 

P s=p 



ei = [1 0 ... 0] r 



and 



a = [1 5(1) .. . a(p) ] T ; (p + 1 elements) 



Equation 4.12 is conceptually solved as follows : 



b = C~ l ei 

b = ^ (4.14) 

CF~ 

Since a(0) = 1, a 2 — and a = <r 2 b. 

Unfortunately C~ l does not exist if C is singular. For this reason, the 
Moore-Penrose pseudoinverse, C + , is used. This may be computed from 
the nonzero eigenvalues, A*, and corresponding eigenvectors, v*, of C as 
follows: 

M 

c + = £ -V.vf (4.15) 

,=i 

where the eigenvalues are ranked in decreasing magnitude Aj > A 2 . . . A p 
and there are m estimated principle components (m < p). The value m 
represents the effective rank of the matrix C- 



ii) Singular' Value Decomposition 

As already stated, we chose not to use the standard normal equation solu- 
tion to the prediction problem because of its known sensitivity to ill con- 
ditioning. The calculations shown in this paper were made using the SVD 
method proposed by Tufts & Kuniarcsan [GS] and subsequently studied by 
Minanii, Kawaka and Minami [39]. The use of the SVD permits the use of 
a regularized version of the Moore Penrose generalized inverse in which the 
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regularization can be controlled by setting the small singular values to zero 
when they are within a threshold, typically a small multiple of machine pre- 
cision - macheps (smallest number represented by a given computer which 
when added to 1 returns something greater than 1). This threshold is a 
function of the SNR, and should be adjusted accordingly. 

The AWVD is implemented by using the singular value decomposition 
(SVD) of the kernel data matrix to compute a reduced rank principle eigen- 
vector approximation to the inverse covariance matrix for use in the AR 
modelling process. Since C can be expressed as the product of the ob- 
servation matrix, I\, and its Ncrmifian transpose, the eigenvectors and 
eigenvalues of C can be determined from a singular value decomposition of 
I<, i.e. 

C=K n K (4.16) 

where I\ is the matrix defined by the elements k(i,j) — I\ t (p 4- i — j), for 
i — 0 . . . N — 1 — p, j = 0 . . . p. Singular value decomposition of I\ gives 

K = ULV" (4.17) 

where U and V are unitary and £ is a diagonal matrix with the elements, 
or(k,k), ( called the singular values of K) ranked in decreasing order such 
that <r(l, 1) > a( 2,2) . . . > cr(p,p). Then the eigenvectors of C correspond 
to the columns of V and the eigenvalues of C are given by the squares of 
the singular values of K. 

iii) The Algorithm: Details of the Implementation 

For a given value of t, eq (10) can be expressed as: — A'a + n = k, with the 
following notation: 

«p(p) 

«i(p) 

n p{?) 

»i (p) 

’ Mp) 

k = : 

. Mp) 
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4(4) .. 


k(p + 1) 
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K is seen to be a Hankel matrix where the elements of the matrix are 
the same along the reverse diagonal and whose first row is kT\ 

In the absence of noise, the estimate a can be given by: 

a = — A' + k. (4.19) 

where K + is the generalized inverse of the matrix K obtained via the 
reduced rank SVD, i.e. is obtained using the inverse of singular values 
above the threshhoid and the corresponding eigenvectors. 

This solution corresponds to a least-squares fit of the data to a model 
using a finite number p of non-harmonically related sinusoids in noise, where 
p is the order of the autoregressive model. The parameter p is estimated 
by either the numerical rank of the data matrix K as computed via the 
SVD, or is estimated by information theoretic methods. The reader is 
asked to compare this result with the observation in [34] that the discrete 
Fourier Transform of a data sequence is equivalent, (up to a normalization 
factor) to a least-squares fit with a model which is composed of a finite 
number ( N , size of the DFT) of harmonically rela ted sinusoids. It is through 
the inclusion of the knowledge that we arc looking for p non harmonically 
non-related sinusoids instead of N harmonically related sinusoids that the 
improved (high-precision and accuracy) performance of this Autoregressive 
Wigner-Ville estimation results. 

Selection of the rank, m, can be made by determining the number of the 
in largest singular values which represent primarily the signal components 
and are usually well separated from the p-m-1 small values which represent 
the noise. In practice, to represent the spectrum of a signal which consists 
of L sinusoids, we need a rank of at least 2 L. Slight overestimation of the 
rank usually does not afTect significantly the resulting spectral estimate. 
The smallest singular values corresponding to the noise (assumed to be 
equal for white noise), determine the value of cr. 

The window length should be chosen such that if is large enough to 
contain the majority of the signal, but bounded by its finite duration. The 
window used is at most the natural effective window — r) of the 

kernel occuring in the WVD process [8]. 
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The cross-section of the Wigner-Ville Distribution at time t is then 
estimated by: 



W(tJ) 



I 1 + EL. b{n)e~ i 2 *J n \ 2 



(4.20) 



where a 2 and d(n) are determined from a kernel, k t , centered about 
time t, where i — n/ j s . 



4.6.2 Implementation of the AWVD 

Conventional WVD analysis as defined in (2) is performed using a software 
package based on the algorithm described in [8]. In the new method, the 
estimated AWVD is calculated as shown in Figure 1 using a peak-finding 
algorithm to ensure accurate representation of the extremely sharp peaks. 
Two important features of the implementation arc: 

1. AWVD performance and speed is improved by performing the mod- 
elling on the real part of the kernel since the imaginary part contains 
redundant information. However the analytic signal must still be used 
to form the kernel, so that low frequency artifacts arc not created [11]. 

2. The analytic signal must be interpolated before forming the kernel. 



4.6.3 Results and discussion 

i) Linear FM 

The signals considered in this section arc FM signals. The linear FM signal 
can be expressed as: 

z(t) = e ;2*Dc<+«< 3 /2) (4.21) 

for which the WVD kernel becomes: 

k,(t,r) = e ' 24 /c+^)r (4.22) 

It is clear that this kernel is composed of one single frequency at any 
particular time, and that a spectral analysis of the windowed kernel will be 
more effective than the spectral analysis of the windowed signal since the 
windowed signal will always contain more than one spectral component. 

The test signal is a 512 sample linear FM signal sampled at 200 II a and 
modulated from a lower frequency of 10 II/. t.o an upper frequency of OOIIz. 
In this implementation, the analytic signal is calculated in the frequency 
domain and is interpolated by zero padding before performing the inverse 
discrete time Fourier Transform [II], The AWVD results arc shown in 
Figure 2 and can be compared with the WVD result in Figure 3. 
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ii) Hyberbolic FM 

The hyperbolic FM signal is expressed as: 

z{t) = e 



(4.23) 



where /o is the initial frequency of the sweep at t = 0 and / c is the 
center frequency. 

Its instantaneous frequency is: 

m = <«•*> 

for which the WVD kernel becomes: 



t) = c 



log{ 



I' + nlM-r/a)! 1 
l+«(l-r/3)| / 



By dividing both numerator and denominator by 2(1 -f at), the expression 
for the log can be expanded as: 



log 



1 + 



ar 



2(1 + al) 



- log 



1 - 



ar 



2(1 + at) 



which can be simplified if we have the condition: <C 1. 

In practice, the signals under analysis are causal and have a finite du- 
ration T. Therefore, the maximum value of the numerator in the above 
inequality is aT. From the expression for /,(£) given in Equation 4.24, we 
can also see that the minimum value of the denominator is 1. The previous 
inequality is therefore valid when aT 2 is satisfied. From the expression 
for fi(t) in (24), when t equals T, then f,(t) equals / 0 + B. This yields the 
relation aT = ^ 2. By introducing the signal center frequency f c , 

this condition can be expressed as: 



B<f c 



(4.25) 



This condition corresponds to narrow band signals. We thus obtain 

K z {L,t) ~ = e ' 2T M0 T (4.26) 

and the same comments as in the linear FM case apply. 
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results 

The AWVD and WVD were applied t.o an increasing hyperbolic FM signal 
and the results are shown in figures 4 and 5. The peaks of the AWVD follow 
the frequency law more closely than the WVD does and the peak width is 
typically reduced by a factor of 5 to 10. Since the Covariance AR modelling 
is not subject to windowing distortion, the multiple sidelobes evident in the 
rectangularly windowed WVD do not appear in the AWVD. Improvement 
in the frequency law tracking of the AWVD over the conventional WVD is 
evident. The autoregressive modelling of the WVD kernel data is better 
than the AR modelling of the signal data because the bilinear form ensures 
that within a window, we have only one frequency (see equation at the 
beginning of this section), whereas with the signal data, we always have 
several frequencies present in the window. 

iii) Multicomponent Signals 

Although it is expected the the AWVD would perform badly for this 
class of signals, because of the cross-terms generated by the bilinear kernel, 
it is nevertheless interesting to present the following example. A multi- 
component signal formed from the sum of two chirp signals is analyzed by 
the AWVD and the WVD. and the results arc shown in fig. 6 and fig. 7. 
This example shows the cross-product “ghost” peaks which alternate be- 
tween positive and negative values in the WVD, but are always positive 
in the AWVD as well as being reduced. The artifacts found in the WVD 
analysis of multicomponent signals are reduced by the use of the AWVD 
because the components of the Kernel which give rise to the artifacts arc 
not well modelled by the Autoregressive model which was chosen to model 
the components of the kernel corresponding to the “autoterms”. 

The robustness of the method to parameter selection and to additive 
noise has been tested [12], For chirp signals analyzed in noise, the results 
are summarized as follows: 

1. For signals with a large SNR, the result is a high accuracy estimate of 
the WVD, which is largely independent, of the model order and rank. 

2. The length of the window appears not to affect the resolution of the 
WVDs. 

3. The rank should be approximately the same as (or greater by 1 than) 
the model order. 

4. The higher the model order and rank the greater the resolution of the 
peaks, but if the model order becomes Loo large, spurious peaks may 



occur. 
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5. The model breaks down below about 5-10 db SNR. This is a reason- 
able result since the method for estimating the WVD peaks is similar 
to Prony’s method for estimating sinusoids in noise - it can only be 
successful if the original signal is dominant, i.e. the SNR is large 
enough. 

4.6.4 Applications of the AWVD 

This section is intended to provide the reader with a few examples of ap- 
plications of this method, which provided the reason for developing this 
AWVD estimator. 

1. Seismic: As the echolocalion search for petroleum resources becomes 
more difficult, the seismic industry must resolve geological structures 
with increasing accuracy. Since the bandwidth of the signals is lim- 
ited by the propagation (attenuation) properties of the earth, mod- 
ern spectral estimation techniques are playing an ever more important 
role in seismic prospecting. The high-accuracy AWVD described here 
provides one method for achieving better identification of fine strata; 
and this can provide enhanced interpretation of seismic sections, re- 
vealing more detail about the thin layers. The modern methods of 
seismic exploration are now based upon the use of linear FM sig- 
nals as source instead of classical explosives. These FM signals are 
emitted by multiple large trucks with vibrating plates in contact with 
the surface of the ground. This signal fits the model used for the 
derivation of the AWVD estimator. Under the assumption of the 
convolution model (weak reflection), then the Wigncr-Ville Distribu- 
tion of the received signal is the convolution (in time) between the 
WVD of the emitted signal and the WVD of the sequence of reflec- 
tors. If the WVD of the sequence of reflectors was known, then using 
the well-known invertability of the WVD, the sequence of reflectors 
can be recovered up to a multiplicative constant. Thus the practical 
problem of Wigner-Villc analysis of the seismic problem is a classical 
inverse scattering problem in which if is desired to deconvolve the 
instrument response function, i.e. the WVD of the transmitted signal 
from the WVD corresponding to the sequence of reflectors. Fourier 
deconvolution procedures are not appropriate in this case, and high 
accuracy Wigner-Villc Analysis using the AWVD may be thought of 
as being equivalent to a numerically robust solution to the classical 
inverse scattering problem. 

2. Sonar/radar problem: The sonar/radar ec.holocation problem is a gen- 
eralization of the seismic problem described above in which the itidi- 
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vidual reflectors may also be in motion so that they both provide a 
delay and a frequency shift of the signal. This problem is normally 
described in terms of the Ambiguity-function of the transmitted wave- 
form which may be viewed as the output of a matched filter matched 
to the signal waveform when that waveform is delayed and frequency- 
shifted. It is well-known [37] that the WVD is equal to the 2DFT of 
the ambiguity function, and therefore the whole echolocation system 
can be described in terms of the WVD [74], In the case where the scat- 
tered signals are modelled as Gaussian random variables, the WVD 
at the output of the receiver may be written as the 2D convolution 
of the WVD of the transmitted signal with the scattering function 
of the reflectors in range and floppier. This, too, is now a 2D scat- 
tering problem where the instrumentation function, ie; the WVD of 
the transmitted signal, approximates a line integral, eg. a projection 
of the scattering distribution when the transmitted signal is a linear 
FM, and thus the equivalent deconvolution problem is anhomorphic 
(there is different accuracy in range and floppier). Also, the deconvo- 
lution is enhanced by high-accuracy techniques which make the line 
integral approximation more closely valid, so that the deconvolution 
can be achieved by tomographic reconstruction techniques [74]. 

3. Fresnel region beamforming: The third example is an example of the 
passive localisation problem where multiple spherical waves are re- 
ceived by a uniformly spaced line-array of sensors such as geophones, 
hydrophones or electromagnetic antennas. It has been shown in the 
array processing literature that when the receiving array is in the 
far field, the spherical waves may be approximated by plane waves. 
Then the distribution of phase along the sensors will be linear and 
different ratios of phase shifts correspond to different angles of ar- 
rival; and thus the problem of estimating the direction of arrival is 
equivalent to the spectral analysis of the time signal corresponding 
to the spatial samples of the signal. This type of signal processing is 
equivalent to decomposing the incoming signal Held into a frequency- 
wavenumber spectrum [28]. When the sources arc in the Fresnel zone 
of the array, then the curvature of the spherical wavefront must be 
properly accounted for. Inspection shows that it is now a linear FM, 
which again corresponds to the signal for which the proposed AWVD 
method will be most effective [13] [66]. Although WVD analysis al- 
lows for the simultaneous estimation of direction of arrival and range, 
a high-accuracy WVD is required to provide practical estimation of 
the range. 
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LIST OF FIGURES AND CAPTIONS 

Figure 1: Time-Frequency Representation of a linear FM chirp signal 
obtained by using an SVD-based Autoregressive Wigner-Villc Distribution 
(AWVD). 

Figure 2: Time-Frequency Representation of a linear FM chirp signal 
obtained by using the conventional Wigner-Ville Distribution (WVD). 

Figure 3: Time-Frequency Representation of a hyperbolic FM chirp 
signal obtained by using an SVD based Autoregressive Wigner-Ville Dis- 
tribution (AWVD). 

Figure 4: Time-Frequency Representation of a hyperbolic FM chirp sig- 
nal obtained by using the conventional Wigner-Ville Distribution (WVD). 

Figure 5: Time-Frequency Representation of two cross-linear FM chirp 
signals obtained by using an SVD based Autoregressive Wigner-Ville Dis- 
tribution. 

Figure 6: Time-Frequency Representation of two cross-linear FM chirp 
signals obtained by using the conventional Wigner-Ville Distribution (WVD) 



4.7 Radar Imaging 

A new formulation of tomographic imaging for rotating targets is presented 
in terms of the Wigner-Ville distribution. This new formulation both sim- 
plifies the signal processing required for image generation and is applicable 
in both mono-static and bi-static situations. Also, the new formulation will 
allow the use of parametric modifications of the WVD, including the AWVD 
described above, to improve the resolution attainable in radar imaging. 

The target is assumed to have multiple scaltcrer within each range- 
Dopplcr resolution cell so that the received signal will be a sample function 
of a Gaussian random variable. When the transmitted waveform is a lin- 
ear frequency modulated chirp then the expected value of the output of 
a matched filter receiver is the two-dimensional convolution of the rangc- 
Doppler scattering function of the target with the aulo-ambiguity surface of 
the transmitted waveform. Under the additional assumption that the trans- 
mitted signal has a Gaussian or rectangular envelope then the receiver can 
be simplified to a short-time Fourier transform receiver using only a prion 
information about the envelope and the chirp rate of the transmitter. This 
simplification eliminates the need for detailed phase information at the re- 
ceiver and thus bi-static operation is feasible. Since most surfaces are rough 
at optical frequencies this form of. signal processing should be partic.ullarly 
useful with optical sensors such as coherent laser radars (LIDARsj. 
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4.7.1 Introduction to Range-Doppler Imaging 

Range-Doppler imaging (scatterer distribution mapping) is a general class 
of inverse problem which arises in remote sensing applications in either 
radar or sonar whenever the sensor is in relative motion witli respect to the 
object or objects under study. In principle the sensors can radiate either 
electromagnetic (radar) or acoustic (sonar) signals. The radiated signals 
are usually in the form of a modulation on a carrier wave. The ratio of 
modulation bandwidth to carrier frequency is called the fractional band- 
width. We can distinguish between two cases, depending on the fractional 
bandwidth: 

1. Small Fractional Bandwidth: The fractional bandwidth is usually small 
in the case of radar and thus motion induced Doppler can be treated 
as a single frequency shift of the entire modulation spectrum. 

2. Large Fractional Bandwidth: The fractional bandwidth is often large 
in the case of sonar and thus motion induced Doppler should be 
treated as a change in time scale of the modulation signal or equiva- 
lently the Doppler shift caused by the motion is proportional to the 
frequency of the components of the transmitted spectrum. 

For simplicity we will consider only the small fractional bandwidth case and 
will use radar terminology and examples. However, the signal processing 
concepts presented in this section can be extended to the large fractional 
bandwidth case for use in sonar signal processing. 

The problem which will be the subject of this paper is the mapping 
of rigid objects which are in motion with respect to a stationary sensor. 
Although rigid objects are the primary targets of interest in most remote 
sensing applications non rigid scatter distributions such as a flight of birds 
or a school of fish may be mapped if the overall shape of the swarm does 
not change too much during the mapping operation. By specalizing to 
rigid objects, it is possible in many applications to incorporate a priori or 
in some cases a posteriori information about the target and thus convert 
the measured scatterer distribution to a range and cross-range image of 
the object. Two types of scattering, deterministic and stochastic, have 
been considered in the literature. This paper will assume the stochastic 
model of scattering which is usually valid when the target is rough with 
respect to the wavelength of the illuminations. Thus the inverse problem 
is to estimate the range-Doppler scattering function of the target which 
represents the average power reflected from a resolution cell of the target. 
If the target is essentially two-dimensional, (i.e. a rotating disk), when 
viewed along the line of sight (LOS) of the radar than the range-Doppler 
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cells of the scattering function form a one-to-one mapping with range and 
cross-range cells of the object. If the target is essentially three-dimensional, 
(i.e. a rotating sphere), then the mapping is no longer one-to-one and in 
general there may be two or more range and cross-range cells with the same 
range and Doppler (e.g. the rotation axis of the sphere is perpendicular to 
the LOS of the radar). 

Although the analysis of the tomographic imaging method in this paper 
is stated in terms of a monostatic radar, (i.e. transmitter and receiver either 
co-located or closely located), the analysis can be extended to a bi-static 
radar, (transmitter and receiver at different locations). In the bi-static 
case the transformation from range and Doppler to range and cross-range 
is more difficult due to the geometrical terms in the transformation [3]. 
However, the measurement of the range-Doppler scattering function is still 
simple when the analysis is performed using the Wigncr-Vilic distribution 
instead of the ambiguity surface, (i.e. magnitude squared of the ambiguity 
function). As will be subsequently shown, the estimation of the scatter- 
ing function separates into a measurement of the Wigner-Ville distribution 
of the received signal made at the receiver and the utilization of a prori 
information about the transmitted signal. 

Although the concept presented in this paper is based on a stochastic 
model of the reflection process, it still works as described even when the 
target is deterministic and not fluctuating. In this latter case the signal 
processing proposed in this paper may be simplified by reducing the num- 
ber of different chirp signals used and thus may be useful at microwave 
frequencies where longer coherent integration times are required due to the 
reduced Doppler sensitivy of the lower carrier frequency, compared to the 
LIDAR case. Also, since the nominal resolution of the range-Doppler cell 
is proportional to the bandwidth (range resolution) and the durat ion of co- 
herent integration (Doppler resolution) use at microwave frequencies may 
require stepped frequency chirp signals instead of continuous chirp signals 
to satisfy the bandwidth limitations of available components [7 1 J. 

4.7.2 WVD Background for Echography 

In order to understand the role of the Wigner-Ville distribution in radar/sonar 
signal processing, it is necessary to introduce a model of the received signal 
which is sufficiently realistic to be practical yet sufficiently simple so that it 
can be easily understood. Therefore, only scalar signals will be considered. 
This assumption is true for acoustic waves propagating in a liquid medium 
and is a simplifying assumption for electromagnetic waves. 

The received signal will be modeled as the linear superposition of sealed 
versions of the transmitted signal delayed by the round trip time and mod- 
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ified by Doppler. Since the transmitted signal has only limited range and 
Doppler resolution then the received signal may be simplified by assuming 
that the scale factor is the product of a constant that depends on the range 
and a complex random variable which acounts for interference between 
unresolved target scatterers and propagation multipaths. A reasonable as- 
sumption, when the individual contributions are sufficiently numerous and 
approximately equal in magnitude with uniformly distributed phase, is that 
the random variable can be modeled as a zero mean Gaussian with variance 
cr(t, f) which depends on the range and Doppler. This function is tradition- 
ally called the scattering function and this target model has been described 
by VanTrees [69]. 

The scattering function will be assumed to be uncorrelated in frequency 
due to orthogonality of the Doppler shifts for narrow band signals. However, 
different delays may be correlated depending on the characteristics of the 
target. In order to be able to include both deterministic as well as stochastic 
variation from pulse to pulse the random variable will be modeled with an 
exponential time correlation which can be generated by a recursive filter. 
Conventional processing only works well when the correlation is high, (i.e. 
the target scattering is essentially deterministic). Speckle noise is associated 
with single look SAR images. 

4.7.3 Comparision of Radar Imaging Techniques 

Four different analyses of the imaging of radar targets have been proposed 
in the literature. 

In the first and simplest the received signal from a spatially distributed 
target is illuminated at successive frequencies as a function of angle. If 
the scattering is assumed to be deterministic then the two-dimensional 
Fourier transform of the received signal gives the range and cross-range 
components of the scattering distribution provided the target is in the far 
field of the illumination and the variation in angle is sufficiently small that 
a rectangular approximation of a polar grid is valid. This type of signal 
processing is often used in radar ranges where the angular orientation can 
be measured and the signal to noise ratio is sufficiently large that matched 
filter processing is not required. 

In the second method matched filtering is used with either continuous or 
stepped frequency chirp signals. Succcsivc range compressed signals from 
a moving target arc aligned in time in order to compensate for the varia- 
tion in range between succesive measurements and then a discrete Fourier 
transform is computed on the complex matched filter outputs for succssive 
measurements. This Fourier transform computes the Doppler frequency as 
a function, of range and is related to the cross-range of a rigid object if the 
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effective azimuth rotation rate relative to the line of sight (LOS) is known. 
This method is effective for deterministic signal returns, such as those from 
metallic aircraft illuminated by microwave frequencies. However, when the 
target is rough relative to the illuminating wavelength then the scattering 
is stochastic and this type of processing is not appropriate. 

In the third method the cross-ambiguity function between the trans- 
mitted and the received signal is computed. For deterministic scattering 
distributions this integral is proportional to the two-dimensional Fourier 
transform of the product of a complex reflection coeficient with the ambi- 
guity function of the transmitted waveform-buL where the proportionality 
factor is complex and range dependant. This method of range- Doppler 
imaging has been analyzed by Feig and Grunbaum and leads to a new form 
of tomographic imaging [22]. 

In the fourth method, the magnitude squared of the cross-ambiguity 
function between the transmitted and the received signal is computed. For 
stochastic scatterers the expected value of the magnitude squared of the 
cross-ambiguity function is proportional to the two-dimensional convolu- 
tion of the range- Doppler scattering function with the magnitude squared 
of the auto-ambiguity function of the transmitted signal. Thus for rotation 
relative to the LOS the effective Doppler of the target becomes coupled to 
the range of the target through the ambiguity function’s magnitude squared 
response. The latter is sometimes referred to as the ambiguity surface. For 
narrowband signals of unit energy the ambiguity surface is everywhere pos- 
itive and has unit volume. Thus it may be considered to be a probability 
density function which distributes the energy from the scattcrcr in both 
range and Doppler acording to the particular wave form used by the trans- 
mitter. This method of imaging is analagous to the medical method of 
time-of-flight positron cminision tomography (TO FP FT) when finite dura- 
tion chirp signals arc used and is analogous to the medical method of com- 
puterized axial tomography (CAT) in the limit of infinite duration chirp 
signals. The analysis of this method of signal processing will be considered 
next. 

4.7.4 Tomographic Radar Imaging using Ambiguity 
Surfaces 

In 1984, Bernfeld proposed the use of a computerized axial tomography 
(CAT) analogy to process linear fin synthetic aperture (SAlt) images [5]. 
Instead of transmitting repeatedly the same linear fm signal as the radar 
moves along its path, Bernfield suggested that chirps of different slopes 
be transmitted. For a stochastic Gaussian assumption of the scattering 
process the output of a matched filter receiver may be modeled as the two- 
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dimensional convolution of the range-Doppler scattering function with the 
auto-ambiguity surface of the transmitted signal. In the limit of an infinite 
duration linear fm signal and no additive noise the output of the receiver 
may be considered to be equivalent to a sequence of line integrals through 
the scattering function where the orientation of the particular line integral 
relative to the scattering distribution is determined by the chirp rate of 
the transmitted signal. The ambiguity surface of an infinite duration linear 
fm signal may be considered to be a delta function with an orientation in 
range and Doppler determined by the chirp rate of the signal. Therefore, the 
convolution of this ambiguity function with the scattering distribution may 
be considered to be a sequence of tomographic projections of the scattering 
distribution. 

In most sonar and radar applications, infinite duration signals are not 
practical.Thereforc, the tomographic signal processing technique proposed 
by Bernfield should be modified to accomodate finite duration signals. For 
a finite duration linear fm signal the ambiguity surface is only an approxi- 
mation to a line integral. In 1985 Snyder and Whitehouse proposed a sig- 
nal processing modification based on the medical method of time-of-flight 
positron emission tomography (TOFPET) [58], 

When a linear fm signal with a Gaussian envelope is used as the trans- 
mitted signal in an echo ranging system then the auto-ambiguity surface 
of the transmitted signal is an asymmetric two-dimensional Gaussian func- 
tion. The inverse problem which results when a number of dilTcrent orien- 
tations of the two- dimensional Gaussian are convolved with an unknown 
but bounded two-dimensional distribution has been studied in the context 
of reconstructing medical radio nuclcide activity distributions [57]. One 
of the techniques of solving this inverse problem is called the “confidence 
weighting” algorithm and is a generalization of the CAT technique of back 
projection and superposition followed by two-dimensional filtering. 

A phcnomological description of confidence weighting can be easily con- 
structed to motivate the confidence weighting tomographic algorithm. In 
the conventional CAT algorithm there is no information in the projection 
or line integral about the location of the contributions to the line integral. 
Therefore, the value of the line integral is uniformly distributed aiong the 
path of the line integral, (i.e. back projection). In the case of TOFPET or 
the proposed tomographic radar imaging algorithm the pointspread func- 
tion of the tomograph or the ambiguity surface provides a jniori information 
about the contributions to the receiver output. Since the normalized am- 
biguity surface has unit volume and is positive, if may be interpreted as a 
probability distribution and thus confidence weighting is simply a method 
analogous to back projection for incorporating the a priori knowledge. 
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4.7.5 WVD Interpretation of Imaging 

It is well known that the Wigner-Ville distribution (WVD) and the ambigu- 
ity function of Sussman arc related by a two-dimensional Fourier transform 
[67]. Thus it is not surprising that radar or sonar detection and estima- 
tion can be reinterpreted in terms of the WVD. However, it is not without 
some difficulty as will be presented in the following section. But the dif- 
ficulties are minor compared to the insight provided by reformulating the 
signal processing equations. In particular, such an interpretation of the 
signal processing for bi-static operation of the radar or sonar is easier to 
understand, since the effect of receiver mismatch may be predicted and 
in the case of linear frequency modulated (LFM) signals simplified signal 
processing either in terms of the short time Fourier transform (ST FT) or 
STRETCH processing is possible. 

The WVD has many interesting properties [37]. However, only a few 
will be needed for our purposes. These arc: 

• Realness 

• The shift properties 

• Moyal’s formula 

• The marginal properties 

• Relationship to the ambiguity function 

• Relationship to the ambiguity surface 

Unfortunately, there are some additional properties which would be desir- 
able in a time-frequency distribution, but which arc attainable only when 
other desirable prorerties arc sacraficed. These arc: 

• Linearity; Must give up realness or resolution. 

• Positivity; Must give up resolution, the marginal properties, and the 
general applicability of Moyal’s formula. 

Tomographic signal reconstruction techniques restores linearity and as 
well as positivity. High resolution Wigner-Ville analysis shows promise 
in restoring the resolution [12]. However, Moyal’s formula is no longer 
applicable, and therefore there is a small loss of detectability which is not 
usually important for imaging applications [23]. 

The tomographic processing algorithm using the Wigner-Ville distribu- 
tion can be simply deduced from the corresponding tomographic algorithm 
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using the ambiguity surface. The expected value of the cross-ambiguity 
surface of the matched filter receiver is the two-dimensional convolution of 
the scattering function with the auto-ambiguity function of the transmit- 
ted waveform [691. Using the relationship between the ambiguity surface 
and the Wigner-Ville distribution the cross- ambiguity surface is rewritten 
as the two-dimensional convolution of the Wigner-Ville distribution of the 
received signal with the Wigner-Ville distribution of the transmitted signal. 
In a simillar manner the auto-ambiguity surface of the transmitted signal 
is the two-dimensional convolution of the Wigner-Ville distribution with 
itself. If the two-dimensional Fourier transform (i.e. ambiguity function) of 
the WVD of the transmitted signal is positive then the an equivalent for- 
mulation is that the expected value of the Wigner-Ville distribution of the 
received signal is the two-dimensional convolution of the target scattering 
function with the Wigner- Villc distribution of the transmitted signal. 

LASER detection and ranging (LIDA It) technology has significant reso- 
lution advantages over conventional microwave radar techniques. The tech- 
nology has received considerable attention as a means of imaging space- 
borne objects using inverse synthetic aperture radar (ISAlt) processing 
techniques using range resolved tomographic reconstruction. These ISAR 
techniques, however, become ineffective when the object is only viewed for 
a small fraction of a rotation, a common occurrence in strategic defense 
applications. Furthermore, the formidable signal processing requirements 
inherent in a ISAR imaging application require innovative techniques to 
minimize the size, weight, and power of the signal processor, particularly for 
spaceborne systems. Recently, Whitehousc and Boashash have suggested 
that a LIDAR imaging processor which uses tomographic reconstruction 
techniques to form an image may overcome the limited viewing angle prob- 
lem associated with range resolved tomographic ISAR processing [75]. In 
this new technique, the chirp rate of the LIDAR transmitted pulse is varied 
during the illumination period. For these signals, the Wigner-Ville distri- 
bution, which simultaneously measures the signal time-delay and Doppler 
shift, is known to be an approximation to a knife-edge whose orientation 
in the time-delay, Doppler plane is a function of the chirp rate. Thus the 
output of a Wigner-Ville receiver which is the convolution of this rotated 
function with the received LIDAR signal and corresponds to calculating a 
projection of the object’s scattering distribution. Time of flight positron 
emmision tomography (TOFPET) reconstruction techniques can then be 
employed to reconstruct an image of the object by confidence weighting, 
(i.e. convolving the received Wigner-Ville distribution with the Wigner- 
Ville distribution of the illuminating signal), and filtering the superposition 
of the confidence weighted signals. 
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4.8 Conclusions 

The direction finding algorithms we have considered impose a substantial 
computational load for their real-time implementation. Most of this com- 
putational load consists of matrix operations: 

• matrix-vector multiplication 

• linear equation solution 

• orthogonal triangularization 

• triangular backsoive 

• Ilermitian symmetric eigensystem solution 

• singular value decomposition 

• Ilermitian symmetric generalized eigensystem solution 

These matrix operations also comprise a substantial part of the computa- 
tional load for several useful methods for data compression, spectral anal- 
ysis, and pattern recognition. 

4.8.1 Research Issues in Parameter Estimation 

Two of most promising parameter estimation techniques for further study 
for application to direction finding are Cadzow’s method and the Maximum 
Likelihood method. MUSIC and MVDR should be viewed as baseline algo- 
rithms for comparison with the newer methods. Both Cadzow’s eigenvector 
fitting method and the Maximum Likelihood Method permit high resolu- 
tion direction finding in the presence of fully coherent multipath, and both 
are applicable to very general array geometries. However, both require the 
solution of a computationally expensive nonlinear least squares problem. 
At present, it is not only difficult to specify a preferred method for solv- 
ing this nonlinear least squares problem, but even to specify the number 
of arithmetic operations needed for an acceptable solution. In order to 
make Cadzow’s method and MLM fcasable for real-time implementation 
and examine their robustness, further study is needed for: 

• Convergence of Nonlinear Least Squares Algorithms 

• Parallelization of Nonlinear Least Squares Algorithms 

• NLS solution accuracy needed for Cadzow’s Method and MLM 
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• Effect of element gain, phase, and position errors on Cadzow’s Method 
and MLM 

For the special case of radar imaging via estimation of the target scatter- 
ing function, time-frequency analysis via the application of signal subspace 
techniques to the YVigner-Ville distribution leads to problems analogous to 
spectral estimation [76]. For such problems, the spectral analysis of the 
WVD’s bilinear kernel is analagous to direction finding using a uniformly 
spaced line array, and ESPRIT and its variants are strong candidate meth- 
ods. These techniques will require parallel implementations for: 

• The Nonsymmctric Generalized Eigenvalue Problem 

• Total beast Squares 

• The Orthogonal Procrustes Problem 

4.8.2 Parallel Computation 

Surveys of parallel algorithms and architectures for many of the matrix 
computations discussed in this paper are available [62], [64] . Although 
nonlinear least squares algorithms use the preceding operations, the result- 
ing iterative algorithms are not as well understood as the classical matrix 
operations. Additional research is also needed in the area of parallel algo- 
rithms for the nonsymmctric eigenvalue problem, since current techniques 
using parallel Jacobi-like methods for computing the Scluir decomposition 
can have convergence problems for non-normal matrices. [20] 
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Figure 1: Time-Frequency Representation of a linear FM chirp signal 
obtained by using an SVD-bascd Autoregressive Wigncr-Villc Distri- 
bution (AWVD). 
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Figure 2: Time-Frequency Representation of a linear FM chirp signal 
obtained by using the conventional Wigncr-Villc Distribution (WVD). 
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Figure 3: Time-Frequency Representation of a hyperbolic FM chirp 
signal obtained by using an SVD based Autoregressive Wigner-Ville 
Distribution (AWVD). 
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Figure 4: Time-Frequency Representation of a hyperbolic FM chirp 
signal obtained by using the conventional Wigner-Villc Distribution 
(WVD). 
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Figure 5: Time-Frequency Representation of two cross-linear FM 
chirp signals obtained by using an SVD based Autoregressive Wigner-Ville 
Distribution. 
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Figure 6: Time-Frequency Representation of two cross-linear FM 
chirp signals obtained by using the conventional Wigner-Ville Dis- 
tribution (WVD). 




